Load the Libraries

library(sf)
## Warning: package 'sf' was built under R version 4.1.3
## Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 7.2.1; sf_use_s2() is TRUE
library(sp)
## Warning: package 'sp' was built under R version 4.1.3
library(gstat)
## Warning: package 'gstat' was built under R version 4.1.3
library(raster)
## Warning: package 'raster' was built under R version 4.1.3
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.1.3
## -- Attaching packages --------------------------------------- tidyverse 1.3.2 --
## v ggplot2 3.3.6     v purrr   0.3.4
## v tibble  3.1.7     v dplyr   1.0.9
## v tidyr   1.2.0     v stringr 1.4.0
## v readr   2.1.2     v forcats 0.5.1
## Warning: package 'ggplot2' was built under R version 4.1.3
## Warning: package 'tibble' was built under R version 4.1.3
## Warning: package 'tidyr' was built under R version 4.1.3
## Warning: package 'readr' was built under R version 4.1.3
## Warning: package 'dplyr' was built under R version 4.1.3
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x tidyr::extract() masks raster::extract()
## x dplyr::filter()  masks stats::filter()
## x dplyr::lag()     masks stats::lag()
## x dplyr::select()  masks raster::select()
library(spacetime)
## Warning: package 'spacetime' was built under R version 4.1.3
library(stringr)
library(readr)
library(leaflet)
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine
library(gifski)
## Warning: package 'gifski' was built under R version 4.1.3
library(nngeo)
## Warning: package 'nngeo' was built under R version 4.1.3
# to make the spherical geometry errors go away
sf::sf_use_s2(FALSE)
## Spherical geometry (s2) switched off

Import Data

Air quality & sensor data

# AIR QUALITY
## corrected air quality
AQ_df = read.csv("../intermediary_outputs/corrected_AQ_data.csv")

## sensor info
sensors = read.csv("../intermediary_outputs/sensor_data_full.csv")

Geometries

# BOUNDARIES
# Boulder boundaries
BO_CO = st_read("../GIS_inputs_destruction_fireboundary/Boulder_county_munis/Municipalities.shp")
## Reading layer `Municipalities' from data source 
##   `C:\Users\zclem\Documents\GitHub\marshall_fires\GIS_inputs_destruction_fireboundary\Boulder_county_munis\Municipalities.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 10 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -105.525 ymin: 39.91405 xmax: -105.0528 ymax: 40.23771
## Geodetic CRS:  WGS 84
fire_counties = BO_CO %>%
  dplyr::filter(ZONEDESC == "Louisville" |
                  ZONEDESC == "Superior" |
                  ZONEDESC == "Broomfield" |
                  ZONEDESC == "Lafayette" |
                  ZONEDESC == "Boulder")
# Boulder Boundary
boulder_precincts = st_read("../GIS_inputs_destruction_fireboundary/Unincorporated_Boulder/Unincorporated_Boulder.shp")
## Reading layer `Unincorporated_Boulder' from data source 
##   `C:\Users\zclem\Documents\GitHub\marshall_fires\GIS_inputs_destruction_fireboundary\Unincorporated_Boulder\Unincorporated_Boulder.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 129 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -105.4297 ymin: 39.91371 xmax: -105.0528 ymax: 40.21256
## Geodetic CRS:  WGS 84
# Broomfield Boundary
broomfield_precincts = st_read("../GIS_inputs_destruction_fireboundary/Broomfield_Precincts/Precincts.shp")
## Reading layer `Precincts' from data source 
##   `C:\Users\zclem\Documents\GitHub\marshall_fires\GIS_inputs_destruction_fireboundary\Broomfield_Precincts\Precincts.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 61 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 3093759 ymin: 1202697 xmax: 3150881 ymax: 1259386
## Projected CRS: NAD83(HARN) / Colorado North (ftUS)
# Westminster Boundary
westminster_city = st_read("../GIS_inputs_destruction_fireboundary/Westminster_CityLimits/CityLimits.shp")
## Reading layer `CityLimits' from data source 
##   `C:\Users\zclem\Documents\GitHub\marshall_fires\GIS_inputs_destruction_fireboundary\Westminster_CityLimits\CityLimits.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 34 features and 18 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -105.1654 ymin: 39.81811 xmax: -104.9877 ymax: 39.96854
## Geodetic CRS:  WGS 84
# make all the precincts have the same columns (we don't need a lot of the data)
boulder_precincts <- boulder_precincts %>%
  select(c(OBJECTID, geometry, SHAPEarea, SHAPElen))

broomfield_precincts <- broomfield_precincts %>%
  mutate(SHAPEarea = Shape__Are,
         SHAPElen = Shape__Len,
         OBJECTID = OBJECTID + 1000) %>% # add 1,000 to the broomfield precincts to get rid of overlapping labels
  select(-c("GlobalID", "GIS_ID", "PRECINCT_N", "MAP_COLOR", "USER_COMME", "LAST_UPDAT", "Shape__Are", "Shape__Len"))

westminster_city <- westminster_city %>%
  mutate(SHAPEarea = ShapeSTAre,
         SHAPElen = ShapeSTLen) %>%
  select(c(OBJECTID, geometry, SHAPEarea, SHAPElen))

# save proj string for fire-affected counties to use for other transformations
prg = raster::crs(fire_counties,asText=TRUE)

# all boundaries joined together
all_bounds <- boulder_precincts %>%
  rbind(st_transform(broomfield_precincts, crs=prg)) %>%
  rbind(st_transform(westminster_city, crs=prg))
# remove the holes that existed so we just have the outside boundary of this area
all_bounds <- all_bounds %>%
  st_union() %>%
  st_remove_holes()
## although coordinates are longitude/latitude, st_union assumes that they are planar
# Marshall fire boundary
wfigs_fire = st_read("../GIS_inputs_destruction_fireboundary/WFIGS_-_Wildland_Fire_Perimeters_Full_History/FH_Perimeter.shp")
## Reading layer `FH_Perimeter' from data source 
##   `C:\Users\zclem\Documents\GitHub\marshall_fires\GIS_inputs_destruction_fireboundary\WFIGS_-_Wildland_Fire_Perimeters_Full_History\FH_Perimeter.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 9500 features and 108 fields (with 4 geometries empty)
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -165.8152 ymin: 13.38021 xmax: 144.6906 ymax: 69.05174
## Geodetic CRS:  WGS 84
# filter fire to marshall fire only; update crs
marshall_fire = wfigs_fire %>% filter(poly_Incid == "Marshall") %>% st_transform(crs = st_crs(prg))
# BUILDINGS
# Read in destroyed home and building sites
destroyed_homes = st_read("../GIS_inputs_destruction_fireboundary/output_damage_files/destroyed_homes.shp")
## Reading layer `destroyed_homes' from data source 
##   `C:\Users\zclem\Documents\GitHub\marshall_fires\GIS_inputs_destruction_fireboundary\output_damage_files\destroyed_homes.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1030 features and 25 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -105.23 ymin: 39.93075 xmax: -105.135 ymax: 40.01238
## Geodetic CRS:  WGS 84
damaged_homes = st_read("../GIS_inputs_destruction_fireboundary/output_damage_files/damaged_homes.shp")
## Reading layer `damaged_homes' from data source 
##   `C:\Users\zclem\Documents\GitHub\marshall_fires\GIS_inputs_destruction_fireboundary\output_damage_files\damaged_homes.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 143 features and 25 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -105.2309 ymin: 39.93081 xmax: -105.1439 ymax: 40.01083
## Geodetic CRS:  WGS 84
destroyed_businesses = st_read("../GIS_inputs_destruction_fireboundary/output_damage_files/destroyed_businesses.shp")
## Reading layer `destroyed_businesses' from data source 
##   `C:\Users\zclem\Documents\GitHub\marshall_fires\GIS_inputs_destruction_fireboundary\output_damage_files\destroyed_businesses.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 7 features and 24 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -105.1651 ymin: 39.94129 xmax: -105.1384 ymax: 39.95788
## Geodetic CRS:  WGS 84
damaged_businesses = st_read("../GIS_inputs_destruction_fireboundary/output_damage_files/damaged_businesses.shp")
## Reading layer `damaged_businesses' from data source 
##   `C:\Users\zclem\Documents\GitHub\marshall_fires\GIS_inputs_destruction_fireboundary\output_damage_files\damaged_businesses.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 27 features and 24 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -105.2132 ymin: 39.94129 xmax: -105.1384 ymax: 40.01153
## Geodetic CRS:  WGS 84
# create 'type' column; select only the columns needed; make projection string match
destroyed_homes = destroyed_homes %>%
  mutate(type = "destroyed - residential") %>%
  dplyr::rename(jurisdiction = JURISDICTI) %>%
  dplyr::select(jurisdiction, type, latlong, geometry) %>%
  st_transform(crs = st_crs(prg))
destroyed_businesses = destroyed_businesses %>%
  mutate(type = "destroyed - non-residential") %>%
  dplyr::rename(jurisdiction = Jurisdicti) %>%
  dplyr::select(jurisdiction, type, latlong, geometry) %>%
  st_transform(crs = st_crs(prg))
damaged_homes = damaged_homes %>%
  mutate(type = "damaged - residential") %>%
  dplyr::rename(jurisdiction = JURISDICTI) %>%
  dplyr::select(jurisdiction, type, latlong, geometry) %>%
  st_transform(crs = st_crs(prg))
damaged_business = damaged_businesses %>%
  mutate(type = "damaged - non-residential") %>%
  dplyr::rename(jurisdiction = Jurisdicti) %>%
  dplyr::select(jurisdiction, type, latlong, geometry) %>%
  st_transform(crs = st_crs(prg))

# combine into one df, using the "," indicator keeps decimal points
(destroyed_damaged =rbind(destroyed_homes,destroyed_businesses,damaged_homes,damaged_business) %>% separate(latlong, c("lat","long"), ",", convert = FALSE))
## Simple feature collection with 1207 features and 4 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -105.2309 ymin: 39.93075 xmax: -105.135 ymax: 40.01238
## CRS:           +proj=longlat +datum=WGS84 +no_defs
## First 10 features:
##    jurisdiction                    type               lat                long
## 1    Louisville destroyed - residential         39.977609         -105.163733
## 2    Louisville destroyed - residential          39.97682         -105.163673
## 3    Louisville destroyed - residential          39.97682         -105.163673
## 4    Louisville destroyed - residential          39.97682         -105.163673
## 5    Louisville destroyed - residential           39.9777         -105.163555
## 6    Louisville destroyed - residential 39.97771569708369 -105.16356686407488
## 7    Louisville destroyed - residential        39.9777426 -105.16414479323862
## 8    Louisville destroyed - residential 39.97769007192905 -105.16382127167078
## 9    Louisville destroyed - residential 39.97774709125106 -105.16359059222464
## 10   Louisville destroyed - residential         39.976826         -105.162003
##                      geometry
## 1  POINT (-105.1637 39.97761)
## 2  POINT (-105.1637 39.97682)
## 3  POINT (-105.1637 39.97682)
## 4  POINT (-105.1637 39.97682)
## 5   POINT (-105.1636 39.9777)
## 6  POINT (-105.1636 39.97772)
## 7  POINT (-105.1641 39.97774)
## 8  POINT (-105.1638 39.97769)
## 9  POINT (-105.1636 39.97775)
## 10  POINT (-105.162 39.97683)
# change lat and long to floats
destroyed_damaged$lat = as.double(destroyed_damaged$lat)
destroyed_damaged$long = as.double(destroyed_damaged$long)

Maps

# create "location status" column for mapping]
sensors$loc_status = str_c(sensors$Location, " - ", sensors$Status)
sensors = na.omit(sensors)
# set factors to false
options(stringsAsFactors = FALSE)

# create palette for destroyed and damaged building colors for leaflet map
pal = colorFactor(c("red","orange","magenta","light pink"), domain = unique(destroyed_damaged$type))


## create indoor/outdoor icons
sensorIcon = awesomeIconList(
  "inside" = makeAwesomeIcon(
    icon = "home",
    iconColor = "white",
    markerColor = "blue",
    library = "fa"
  ),
  "outside" = makeAwesomeIcon(
    icon = "tree",
    iconColor = "white",
    markerColor = "green",
    library = "fa"
  ))

sensors
##       X     ID fire_period post_fire_period pre_fire_period
## 1     1   3062  100.000000        100.00000       100.00000
## 2     2   4583  100.000000         99.93082       100.00000
## 3     3   4643  100.000000         99.33126       100.00000
## 4     4  15663  100.000000         99.52727       100.00000
## 5     5  27120   48.333333         99.86740       100.00000
## 6     6  30743   74.000000         73.96770       100.00000
## 7     7  33109   53.333333         99.70022       100.00000
## 8     8  38623    7.333333         95.18621       100.00000
## 10   10  45325   99.666667         88.22783       100.00000
## 11   11  45369   99.666667         99.94235       100.00000
## 13   13  50133  100.000000        100.00000       100.00000
## 14   14  53729  100.000000         99.84434       100.00000
## 16   16  56209  100.000000         94.93255        96.66667
## 17   17  56431   86.666667         99.87893       100.00000
## 18   18  60507  100.000000         99.93188       100.00000
## 19   19  60517  100.000000         99.82705       100.00000
## 20   20  64671  100.000000         99.81552       100.00000
## 21   21  64967   22.333333         99.86740       100.00000
## 22   22  65669   17.666667         99.89046       100.00000
## 23   23  66173  100.000000         99.96541       100.00000
## 24   24  67917   99.333333         99.96541       100.00000
## 25   25  68375    7.333333         99.25055       100.00000
## 26   26  69103    4.000000         99.72328       100.00000
## 28   28  70241  100.000000         99.96541       100.00000
## 29   29  70911  100.000000         94.94408       100.00000
## 31   31  72931  100.000000         99.94811       100.00000
## 34   34  79121    7.333333         92.70725       100.00000
## 35   35  79907   38.000000         99.99423       100.00000
## 36   36  81149  100.000000         99.59068       100.00000
## 37   37  82879  100.000000         99.68292       100.00000
## 38   38  83423  100.000000         98.81817       100.00000
## 39   39  83597   97.333333         99.58492       100.00000
## 40   40  83661  100.000000         95.07091       100.00000
## 41   41  83785    4.000000         99.72328       100.00000
## 42   42  84169  100.000000         99.91929       100.00000
## 43   43  84585   76.000000         98.76780       100.00000
## 44   44  84697    1.333333         96.62170       100.00000
## 45   45  84777   10.666667         63.97674       100.00000
## 46   46  85643  100.000000        100.00000       100.00000
## 47   47  85999   14.000000         97.89000       100.00000
## 48   48  86477   13.666667         98.16672       100.00000
## 49   49  86637   99.666667         99.78669       100.00000
## 50   50  87465  100.000000         99.98847       100.00000
## 51   51  87531  100.000000         99.38314       100.00000
## 52   52  88035  100.000000         99.52150       100.00000
## 53   53  88533  100.000000         99.97694       100.00000
## 55   55  90255    3.666667         99.64257       100.00000
## 56   56  90297   89.000000         90.30900       100.00000
## 57   57  91877  100.000000         99.94811       100.00000
## 58   58  92321  100.000000         99.97117       100.00000
## 59   59  94411  100.000000         88.86199       100.00000
## 60   60  95591  100.000000        100.00000       100.00000
## 61   61  95831   99.000000         87.35155       100.00000
## 62   62  96731  100.000000        100.00000       100.00000
## 63   63  97647  100.000000         99.93658       100.00000
## 64   64 100521  100.000000         99.60798       100.00000
## 65   65 101003    4.000000         99.40044       100.00000
## 66   66 102614  100.000000         99.47538       100.00000
## 67   67 102724  100.000000         99.82128       100.00000
## 69   69 103724    7.333333         98.20708       100.00000
## 70   70 103820   49.000000         46.41992       100.00000
## 71   71 104594   99.666667         99.78093       100.00000
## 72   72 106214   22.333333         99.85587       100.00000
## 73   73 106492   48.666667         98.85276       100.00000
## 74   74 106800  100.000000         99.67139       100.00000
## 75   75 106906   99.666667         99.99423       100.00000
## 76   76 107466   63.000000         99.96541       100.00000
## 78   78 108558   99.666667         99.77516       100.00000
## 80   80 110042  100.000000        100.00000       100.00000
## 81   81 110122   20.666667         99.96541       100.00000
## 82   82 111834   56.666667         99.34855       100.00000
## 83   83 112082   56.666667         99.59068       100.00000
## 84   84 112522  100.000000         97.51528       100.00000
## 85   85 112606   99.666667         97.58446       100.00000
## 86   86 112974  100.000000         99.98847       100.00000
## 87   87 113054  100.000000         99.96541       100.00000
## 88   88 113266  100.000000         99.87317       100.00000
## 93   93 115303   99.666667         99.94811       100.00000
## 94   94 117647  100.000000         99.82705       100.00000
## 95   95 119547   99.666667         99.44656       100.00000
## 97   97 120477  100.000000         99.96541       100.00000
## 99   99 120847   53.333333         99.98812       100.00000
## 101 101 121197  100.000000         99.98656       100.00000
## 103 103 122071   98.333333         99.00842       100.00000
## 106 106 122585    4.000000         96.54675       100.00000
## 115 115 124319   53.333333         99.46962       100.00000
## 200 200 141956   45.000000         64.57397       100.00000
##                                                Status                    Month
## 1                Complete data throughout fire period Before or during 12-2021
## 2                Complete data throughout fire period Before or during 12-2021
## 3                Complete data throughout fire period Before or during 12-2021
## 4                Complete data throughout fire period Before or during 12-2021
## 5         Sensor offline during fire, returned online Before or during 12-2021
## 6   Sensor offline during fire, did not return online Before or during 12-2021
## 7         Sensor offline during fire, returned online Before or during 12-2021
## 8         Sensor offline during fire, returned online Before or during 12-2021
## 10               Complete data throughout fire period Before or during 12-2021
## 11               Complete data throughout fire period Before or during 12-2021
## 13               Complete data throughout fire period Before or during 12-2021
## 14               Complete data throughout fire period Before or during 12-2021
## 16               Complete data throughout fire period Before or during 12-2021
## 17               Complete data throughout fire period Before or during 12-2021
## 18               Complete data throughout fire period Before or during 12-2021
## 19               Complete data throughout fire period Before or during 12-2021
## 20               Complete data throughout fire period Before or during 12-2021
## 21        Sensor offline during fire, returned online Before or during 12-2021
## 22        Sensor offline during fire, returned online Before or during 12-2021
## 23               Complete data throughout fire period Before or during 12-2021
## 24               Complete data throughout fire period Before or during 12-2021
## 25        Sensor offline during fire, returned online Before or during 12-2021
## 26        Sensor offline during fire, returned online Before or during 12-2021
## 28               Complete data throughout fire period Before or during 12-2021
## 29               Complete data throughout fire period Before or during 12-2021
## 31               Complete data throughout fire period Before or during 12-2021
## 34        Sensor offline during fire, returned online Before or during 12-2021
## 35        Sensor offline during fire, returned online Before or during 12-2021
## 36               Complete data throughout fire period Before or during 12-2021
## 37               Complete data throughout fire period Before or during 12-2021
## 38               Complete data throughout fire period Before or during 12-2021
## 39               Complete data throughout fire period Before or during 12-2021
## 40               Complete data throughout fire period Before or during 12-2021
## 41        Sensor offline during fire, returned online Before or during 12-2021
## 42               Complete data throughout fire period Before or during 12-2021
## 43               Complete data throughout fire period Before or during 12-2021
## 44        Sensor offline during fire, returned online Before or during 12-2021
## 45  Sensor offline during fire, did not return online Before or during 12-2021
## 46               Complete data throughout fire period Before or during 12-2021
## 47        Sensor offline during fire, returned online Before or during 12-2021
## 48        Sensor offline during fire, returned online Before or during 12-2021
## 49               Complete data throughout fire period Before or during 12-2021
## 50               Complete data throughout fire period Before or during 12-2021
## 51               Complete data throughout fire period Before or during 12-2021
## 52               Complete data throughout fire period Before or during 12-2021
## 53               Complete data throughout fire period Before or during 12-2021
## 55        Sensor offline during fire, returned online Before or during 12-2021
## 56               Complete data throughout fire period Before or during 12-2021
## 57               Complete data throughout fire period Before or during 12-2021
## 58               Complete data throughout fire period Before or during 12-2021
## 59               Complete data throughout fire period Before or during 12-2021
## 60               Complete data throughout fire period Before or during 12-2021
## 61               Complete data throughout fire period Before or during 12-2021
## 62               Complete data throughout fire period Before or during 12-2021
## 63               Complete data throughout fire period Before or during 12-2021
## 64               Complete data throughout fire period Before or during 12-2021
## 65        Sensor offline during fire, returned online Before or during 12-2021
## 66               Complete data throughout fire period Before or during 12-2021
## 67               Complete data throughout fire period Before or during 12-2021
## 69        Sensor offline during fire, returned online Before or during 12-2021
## 70  Sensor offline during fire, did not return online Before or during 12-2021
## 71               Complete data throughout fire period Before or during 12-2021
## 72        Sensor offline during fire, returned online Before or during 12-2021
## 73        Sensor offline during fire, returned online Before or during 12-2021
## 74               Complete data throughout fire period Before or during 12-2021
## 75               Complete data throughout fire period Before or during 12-2021
## 76        Sensor offline during fire, returned online Before or during 12-2021
## 78               Complete data throughout fire period Before or during 12-2021
## 80               Complete data throughout fire period Before or during 12-2021
## 81        Sensor offline during fire, returned online Before or during 12-2021
## 82        Sensor offline during fire, returned online Before or during 12-2021
## 83        Sensor offline during fire, returned online Before or during 12-2021
## 84               Complete data throughout fire period Before or during 12-2021
## 85               Complete data throughout fire period Before or during 12-2021
## 86               Complete data throughout fire period Before or during 12-2021
## 87               Complete data throughout fire period Before or during 12-2021
## 88               Complete data throughout fire period Before or during 12-2021
## 93               Complete data throughout fire period Before or during 12-2021
## 94               Complete data throughout fire period Before or during 12-2021
## 95               Complete data throughout fire period Before or during 12-2021
## 97               Complete data throughout fire period Before or during 12-2021
## 99        Sensor offline during fire, returned online Before or during 12-2021
## 101              Complete data throughout fire period Before or during 12-2021
## 103              Complete data throughout fire period Before or during 12-2021
## 106       Sensor offline during fire, returned online Before or during 12-2021
## 115       Sensor offline during fire, returned online Before or during 12-2021
## 200 Sensor offline during fire, did not return online Before or during 12-2021
##           Lon      Lat                                        Name Location
## 1   -105.2713 40.04077                                      Kalmia  outside
## 2   -105.1680 40.13818                     Clover Basin Consulting  outside
## 3   -105.2515 39.98943                                       JENSA  outside
## 4   -105.1550 40.10230                         FRCC Niwot location  outside
## 5   -105.2600 39.98123                            Table Mesa Court  outside
## 6   -105.1283 40.13157                                   Creekside  outside
## 7   -105.2622 39.98305                                        Bear  outside
## 8   -105.1749 39.98661                            76th & S Bldr Rd  outside
## 10  -105.1647 40.06149                                Kincross Way  outside
## 11  -105.1649 40.06169                                Kincross Way   inside
## 13  -105.1727 40.10403                                   HRNS-TEST   inside
## 14  -105.0457 39.99003                                Anthem Ranch  outside
## 16  -105.2691 39.97224                             Stony Hill Road  outside
## 17  -105.2821 40.03435                          Grape_and_Broadway  outside
## 18  -105.2694 40.01557                                       Bella  outside
## 19  -105.1199 39.87919                               STANDLEY LAKE  outside
## 20  -105.2709 40.04019                                      K2055i   inside
## 21  -105.2458 39.98447               Majestic Heights, Boulder, CO  outside
## 22  -105.1477 39.97686                                  Louisville  outside
## 23  -105.1214 39.94361             Terracina Broomfield Apartments  outside
## 24  -105.2833 40.00194                          Chautauqua Heights  outside
## 25  -105.1500 39.98119                       Louisville Air Sensor  outside
## 26  -105.3219 40.04352                                   Hawk Lane  outside
## 28  -105.2804 40.02177                                     MKHOME2  outside
## 29  -105.2508 39.97368                              Shanahan Ridge   inside
## 31  -105.2279 39.99787                               AZTEC CENTRAL  outside
## 34  -105.1741 39.98653                                     Kitchen   inside
## 35  -105.2716 40.03448                                 Vista Drive   inside
## 36  -105.1549 39.97095                                       ntsky  outside
## 37  -105.2401 39.98775                            Mango Manor Coop   inside
## 38  -105.2849 40.05431             Shining Mountain Waldorf School  outside
## 39  -105.2563 39.97676                               South Boulder  outside
## 40  -105.1233 40.13444                        Creekside and Sunset  outside
## 41  -105.3218 40.04347                                   Hawk Lane   inside
## 42  -105.2312 40.01628                                    Snowcave   inside
## 43  -105.2308 39.99910                         645 Manhattan Place   inside
## 44  -105.2798 39.93268                            Eldorado Springs  outside
## 45  -105.1660 39.95349                          Superior Town Hall  outside
## 46  -105.3186 40.02700                           Lower Seven Hills  outside
## 47  -105.1554 39.93897                       Superior - North Pool  outside
## 48  -105.1552 39.92155                       Superior - South Pool  outside
## 49  -105.2515 39.97368                                  Viele Lake  outside
## 50  -105.2845 40.06458                                Dakota Ridge  outside
## 51  -105.2394 39.98681                  The Haunted Pickle Mansion   inside
## 52  -105.1928 40.07116      Gunbarrel Green Spotted Horse (inside)   inside
## 53  -105.1489 39.99006                        Purple house indoors   inside
## 55  -105.2551 39.98033                                   Miami Way  outside
## 56  -105.0926 39.96745                                South Pointe  outside
## 57  -105.1451 39.98836                               Eisenhower Dr  outside
## 58  -105.1491 39.99013                                Purple House  outside
## 59  -105.2689 40.03923                                    Parkside  outside
## 60  -105.2893 40.05224                             Wonderland Lake  outside
## 61  -105.1641 40.06173                              Gunbarrel Hill  outside
## 62  -105.2875 40.02247                                 Concord Ave  outside
## 63  -105.0223 39.92575                     Broomfield - Calkins Pl  outside
## 64  -105.1932 40.07457                                  PurpleHome   inside
## 65  -105.3151 40.04732                                       House   inside
## 66  -105.0764 40.04520                                   Ponderosa  outside
## 67  -105.1251 40.02127                          Blue Heron Estates  outside
## 69  -105.1499 39.97538                                        Home   inside
## 70  -105.2813 40.00443                              Casper’s House   inside
## 71  -105.2007 40.07604                            Bean Mountain Ln  outside
## 72  -105.2458 39.98451                            Majestic Heights   inside
## 73  -105.3364 40.07924                             Boulder Heights   inside
## 74  -105.1249 40.02118                          Blue Heron Estates   inside
## 75  -105.2670 39.97523 1333 Wildwood Ct., Boulder, CO 80305 Inside   inside
## 76  -105.2559 39.98311                              Yale Road Home   inside
## 78  -105.1720 40.09753                               CountryCreek2  outside
## 80  -105.2562 40.03099                            Valmont and 29th  outside
## 81  -105.2427 39.98229                                Hanover Ave.   inside
## 82  -105.2451 40.09622                                    Snead Ct  outside
## 83  -105.2549 40.09183                                   Eagle Ct.  outside
## 84  -105.0474 39.82646                         Linda Park Addition  outside
## 85  -105.2670 39.97515        1333 Wildwood Ct., Boulder, CO 80305  outside
## 86  -105.1615 40.13421                              Summerlin_Lane  outside
## 87  -105.2837 40.04738                                  Quince Ave   inside
## 88  -105.2803 40.04895                         Redwood & Riverside  outside
## 93  -105.0864 39.90192                                 Amli Arista  outside
## 94  -105.1197 39.87945                                Stanley Lake   inside
## 95  -105.0439 39.95181                                  Broadlands  outside
## 97  -105.1246 39.99547                              Bohannan house  outside
## 99  -105.2659 39.97925                                        Yard  outside
## 101 -105.2376 40.00849                               Hancock Drive  outside
## 103 -105.1962 40.04964                                    EMSonJay  outside
## 106 -105.2556 39.98821                                       Bates   inside
## 115 -105.2567 39.98218                              Hartford Drive   inside
## 200 -105.1550 39.97581                          Nighthawk Resident   inside
##                                                                                                                                             address
## 1                                                                      2065, Kalmia Avenue, Boulder, Boulder County, Colorado, 80304, United States
## 2                                                                     5017, William Place, Longmont, Boulder County, Colorado, 80503, United States
## 3                                                                     3565, Eastman Avenue, Boulder, Boulder County, Colorado, 80305, United States
## 4                                                              7037, Longview Drive, Boulder, Niwot, Boulder County, Colorado, 80503, United States
## 5                                                                   2799, Table Mesa Court, Boulder, Boulder County, Colorado, 80305, United States
## 6                                                                  2443, Eagleview Circle, Longmont, Boulder County, Colorado, 80504, United States
## 7                                                                        732, Ithaca Drive, Boulder, Boulder County, Colorado, 80305, United States
## 8                                                                 7618, South Boulder Road, Boulder, Boulder County, Colorado, 80303, United States
## 10                                                         8059, Kincross Way, Heatherwood, Boulder, Boulder County, Colorado, 80301, United States
## 11                                                         8059, Kincross Way, Heatherwood, Boulder, Boulder County, Colorado, 80301, United States
## 13                                                     My Mama's Pie, Murray Street, Boulder, Niwot, Boulder County, Colorado, 80544, United States
## 14                                                                   4512, Silver Mountain Loop, Anthem, Broomfield, Colorado, 80032, United States
## 16                                                                   1916, Stony Hill Road, Boulder, Boulder County, Colorado, 80305, United States
## 17                                      Broadway & Grape Ave, Broadway, Washington Village, Boulder, Boulder County, Colorado, 80306, United States
## 18                                                                      1933, Grove Street, Boulder, Boulder County, Colorado, 80302, United States
## 19                                                    10098, Oak Court, Walnut Grove, Westminster, Jefferson County, Colorado, 80021, United States
## 20                                                                     2099, Kalmia Avenue, Boulder, Boulder County, Colorado, 80304, United States
## 21                                                                  605, South 42nd Street, Boulder, Boulder County, Colorado, 80305, United States
## 22                                                                661, Cleveland Avenue, Louisville, Boulder County, Colorado, 80027, United States
## 23                                                     Building 30, 13630, Via Varra, Overlook District, Broomfield, Colorado, 80020, United States
## 24                                  847, Cascade Avenue, Central Boulder - University Hill, Boulder, Boulder County, Colorado, 80302, United States
## 25                                                            599, West Sagebrush Court, Louisville, Boulder County, Colorado, 80027, United States
## 26                                                                                    67, Hawk Lane, Boulder County, Colorado, 80304, United States
## 28                                                   2357, 13th Street, Washington Village, Boulder, Boulder County, Colorado, 80304, United States
## 29                                                                                          Boulder, Boulder County, Colorado, 80305, United States
## 31                                                                        556, Aztec Drive, Boulder, Boulder County, Colorado, 80303, United States
## 34                                                                      South Boulder Road, Boulder, Boulder County, Colorado, 80027, United States
## 35                                                                       1970, Vista Drive, Boulder, Boulder County, Colorado, 80304, United States
## 36                                                                 775, West Lois Court, Louisville, Boulder County, Colorado, 80027, United States
## 37                                                                      4662, Ingram Court, Boulder, Boulder County, Colorado, 80305, United States
## 38                                     Shining Mountain Waldorf School, 999, Violet Avenue, Boulder, Boulder County, Colorado, 80304, United States
## 39                                                                      1264, Ithaca Drive, Boulder, Boulder County, Colorado, 80305, United States
## 40                                                                     1839, Caleta Trail, Longmont, Boulder County, Colorado, 80504, United States
## 41                                                                                    67, Hawk Lane, Boulder County, Colorado, 80304, United States
## 42                                                                   1707, Commerce Street, Boulder, Boulder County, Colorado, 80301, United States
## 43                                                                    645, Manhattan Drive, Boulder, Boulder County, Colorado, 80303, United States
## 44                                          Eldorado Springs Pool, Artesian Drive, Eldorado Springs, Boulder County, Colorado, 80025, United States
## 45                       Superior Municipal Court, East William Street, Original Superior, Superior, Boulder County, Colorado, 80027, United States
## 46                                                              442, Seven Hills Drive, Seven Hills, Boulder County, Colorado, 80304, United States
## 47  Town of Superior North Pool Facility, South Indiana Street, Rock Creek Ranch, Boulder, Superior, Boulder County, Colorado, 80027, United States
## 48                              Superior South Pool, Huron Peak Drive, Rock Creek Ranch II, Boulder, Boulder County, Colorado, 80027, United States
## 49                                                                     3270, Everett Drive, Boulder, Boulder County, Colorado, 80305, United States
## 50                                                                                   4949, Broadway, Boulder County, Colorado, 80304, United States
## 51                                                                      4682, Ingram Court, Boulder, Boulder County, Colorado, 80305, United States
## 52                                                               5310, Spotted Horse Trail, Boulder, Boulder County, Colorado, 80301, United States
## 53                                                               1998, Eisenhower Drive, Louisville, Boulder County, Colorado, 80027, United States
## 55                                                                          950, Miami Way, Boulder, Boulder County, Colorado, 80305, United States
## 56                                                  South Point Drive, South Pointe, Lafayette, Boulder County, Colorado, 80026-2872, United States
## 57                                                                370, Eisenhower Drive, Louisville, Boulder County, Colorado, 80027, United States
## 58                                                               1998, Eisenhower Drive, Louisville, Boulder County, Colorado, 80027, United States
## 59                                                                             22nd Street, Boulder, Boulder County, Colorado, 80304, United States
## 60                                                                       501, Utica Avenue, Boulder, Boulder County, Colorado, 80303, United States
## 61                                                         8067, Kincross Way, Heatherwood, Boulder, Boulder County, Colorado, 80301, United States
## 62                                                       Concord Alley, Washington Village, Boulder, Boulder County, Colorado, 80302, United States
## 63                                                                                  2830, Calkins Place, Broomfield, Colorado, 80020, United States
## 64                                                                     5583, Homestead Way, Boulder, Boulder County, Colorado, 80301, United States
## 65                                                                230, Timber Lane, Pine Brook Hill, Boulder County, Colorado, 80304, United States
## 66                                                              11821, Flatiron Drive, Leyner, Erie, Boulder County, Colorado, 80026, United States
## 67                                                   2432, Ginny Way, Blue Heron Estates, Lafayette, Boulder County, Colorado, 80026, United States
## 69                                                                   681, Spruce Circle, Louisville, Boulder County, Colorado, 80027, United States
## 70                                   941, Lincoln Place, Central Boulder - University Hill, Boulder, Boulder County, Colorado, 80302, United States
## 71                                                                6655, Bean Mountain Lane, Boulder, Boulder County, Colorado, 80301, United States
## 72                                                                  605, South 42nd Street, Boulder, Boulder County, Colorado, 80305, United States
## 73                                                                 253, Deer Trail Road, Lazy Acres, Boulder County, Colorado, 80455, United States
## 74                                                   2430, Ginny Way, Blue Heron Estates, Lafayette, Boulder County, Colorado, 80026, United States
## 75                                                                    1313, Wildwood Court, Boulder, Boulder County, Colorado, 80305, United States
## 76                                                                          720, Yale Road, Boulder, Boulder County, Colorado, 80305, United States
## 78                                                        7777, Country Creek Drive, Boulder, Niwot, Boulder County, Colorado, 80503, United States
## 80                                                                       3185, 29th Street, Boulder, Boulder County, Colorado, 80301, United States
## 81                                                                  770, South 45th Street, Boulder, Boulder County, Colorado, 80305, United States
## 82                                                                                6732, Snead Court, Boulder County, Colorado, 80503, United States
## 83                                                                                6430, Eagle Court, Boulder County, Colorado, 80503, United States
## 84                                                                    7157, Winona Court, Westminster, Adams County, Colorado, 80030, United States
## 85                                                                    1313, Wildwood Court, Boulder, Boulder County, Colorado, 80305, United States
## 86                                                                                         Longmont, Boulder County, Colorado, 80503, United States
## 87                                                                     1052, Quince Avenue, Boulder, Boulder County, Colorado, 80304, United States
## 88                                                                    1355, Redwood Avenue, Boulder, Boulder County, Colorado, 80304, United States
## 93                                                                          7971, Uptown Avenue, Arista, Broomfield, Colorado, 80021, United States
## 94                                                                 10081, Oak Street, Westminster, Jefferson County, Colorado, 80021, United States
## 95                                                                      4342, Augusta Drive, Broadlands, Broomfield, Colorado, 80023, United States
## 97                           384 Driftwood Circle, 384, Driftwood Circle, Waneka Landing, Lafayette, Boulder County, Colorado, 80026, United States
## 99                                                      Magic Mansion, 2267, Holyoke Drive, Boulder, Boulder County, Colorado, 80305, United States
## 101                                                                    3878, Hancock Drive, Boulder, Boulder County, Colorado, 80303, United States
## 103                                                                               Jay Road, Boulder, Boulder County, Colorado, 80301, United States
## 106                                                                      350, Bates Avenue, Boulder, Boulder County, Colorado, 80305, United States
## 115                                                                    850, Hartford Drive, Boulder, Boulder County, Colorado, 80305, United States
## 200                                                               759, Nighthawk Circle, Louisville, Boulder County, Colorado, 80027, United States
##     zip_code city_update
## 1      80304     Boulder
## 2      80503    Longmont
## 3      80305     Boulder
## 4      80503       Niwot
## 5      80305     Boulder
## 6      80504    Longmont
## 7      80305     Boulder
## 8      80303     Boulder
## 10     80301     Boulder
## 11     80301     Boulder
## 13     80544       Niwot
## 14     80032      States
## 16     80305     Boulder
## 17     80306     Boulder
## 18     80302     Boulder
## 19     80021      States
## 20     80304     Boulder
## 21     80305     Boulder
## 22     80027  Louisville
## 23     80020      States
## 24     80302     Boulder
## 25     80027  Louisville
## 26     80304        Lane
## 28     80304     Boulder
## 29     80305     Boulder
## 31     80303     Boulder
## 34     80027     Boulder
## 35     80304     Boulder
## 36     80027  Louisville
## 37     80305     Boulder
## 38     80304     Boulder
## 39     80305     Boulder
## 40     80504    Longmont
## 41     80304        Lane
## 42     80301     Boulder
## 43     80303     Boulder
## 44     80025     Springs
## 45     80027    Superior
## 46     80304       Hills
## 47     80027    Superior
## 48     80027     Boulder
## 49     80305     Boulder
## 50     80304     Boulder
## 51     80305     Boulder
## 52     80301     Boulder
## 53     80027  Louisville
## 55     80305     Boulder
## 56      2872   Lafayette
## 57     80027  Louisville
## 58     80027  Louisville
## 59     80304     Boulder
## 60     80303     Boulder
## 61     80301     Boulder
## 62     80302     Boulder
## 63     80020      States
## 64     80301     Boulder
## 65     80304        Hill
## 66     80026        Erie
## 67     80026   Lafayette
## 69     80027  Louisville
## 70     80302     Boulder
## 71     80301     Boulder
## 72     80305     Boulder
## 73     80455       Acres
## 74     80026   Lafayette
## 75     80305     Boulder
## 76     80305     Boulder
## 78     80503       Niwot
## 80     80301     Boulder
## 81     80305     Boulder
## 82     80503       Court
## 83     80503       Court
## 84     80030      States
## 85     80305     Boulder
## 86     80503    Longmont
## 87     80304     Boulder
## 88     80304     Boulder
## 93     80021      States
## 94     80021      States
## 95     80023      States
## 97     80026   Lafayette
## 99     80305     Boulder
## 101    80303     Boulder
## 103    80301     Boulder
## 106    80305     Boulder
## 115    80305     Boulder
## 200    80027  Louisville
##                                                      loc_status
## 1                outside - Complete data throughout fire period
## 2                outside - Complete data throughout fire period
## 3                outside - Complete data throughout fire period
## 4                outside - Complete data throughout fire period
## 5         outside - Sensor offline during fire, returned online
## 6   outside - Sensor offline during fire, did not return online
## 7         outside - Sensor offline during fire, returned online
## 8         outside - Sensor offline during fire, returned online
## 10               outside - Complete data throughout fire period
## 11                inside - Complete data throughout fire period
## 13                inside - Complete data throughout fire period
## 14               outside - Complete data throughout fire period
## 16               outside - Complete data throughout fire period
## 17               outside - Complete data throughout fire period
## 18               outside - Complete data throughout fire period
## 19               outside - Complete data throughout fire period
## 20                inside - Complete data throughout fire period
## 21        outside - Sensor offline during fire, returned online
## 22        outside - Sensor offline during fire, returned online
## 23               outside - Complete data throughout fire period
## 24               outside - Complete data throughout fire period
## 25        outside - Sensor offline during fire, returned online
## 26        outside - Sensor offline during fire, returned online
## 28               outside - Complete data throughout fire period
## 29                inside - Complete data throughout fire period
## 31               outside - Complete data throughout fire period
## 34         inside - Sensor offline during fire, returned online
## 35         inside - Sensor offline during fire, returned online
## 36               outside - Complete data throughout fire period
## 37                inside - Complete data throughout fire period
## 38               outside - Complete data throughout fire period
## 39               outside - Complete data throughout fire period
## 40               outside - Complete data throughout fire period
## 41         inside - Sensor offline during fire, returned online
## 42                inside - Complete data throughout fire period
## 43                inside - Complete data throughout fire period
## 44        outside - Sensor offline during fire, returned online
## 45  outside - Sensor offline during fire, did not return online
## 46               outside - Complete data throughout fire period
## 47        outside - Sensor offline during fire, returned online
## 48        outside - Sensor offline during fire, returned online
## 49               outside - Complete data throughout fire period
## 50               outside - Complete data throughout fire period
## 51                inside - Complete data throughout fire period
## 52                inside - Complete data throughout fire period
## 53                inside - Complete data throughout fire period
## 55        outside - Sensor offline during fire, returned online
## 56               outside - Complete data throughout fire period
## 57               outside - Complete data throughout fire period
## 58               outside - Complete data throughout fire period
## 59               outside - Complete data throughout fire period
## 60               outside - Complete data throughout fire period
## 61               outside - Complete data throughout fire period
## 62               outside - Complete data throughout fire period
## 63               outside - Complete data throughout fire period
## 64                inside - Complete data throughout fire period
## 65         inside - Sensor offline during fire, returned online
## 66               outside - Complete data throughout fire period
## 67               outside - Complete data throughout fire period
## 69         inside - Sensor offline during fire, returned online
## 70   inside - Sensor offline during fire, did not return online
## 71               outside - Complete data throughout fire period
## 72         inside - Sensor offline during fire, returned online
## 73         inside - Sensor offline during fire, returned online
## 74                inside - Complete data throughout fire period
## 75                inside - Complete data throughout fire period
## 76         inside - Sensor offline during fire, returned online
## 78               outside - Complete data throughout fire period
## 80               outside - Complete data throughout fire period
## 81         inside - Sensor offline during fire, returned online
## 82        outside - Sensor offline during fire, returned online
## 83        outside - Sensor offline during fire, returned online
## 84               outside - Complete data throughout fire period
## 85               outside - Complete data throughout fire period
## 86               outside - Complete data throughout fire period
## 87                inside - Complete data throughout fire period
## 88               outside - Complete data throughout fire period
## 93               outside - Complete data throughout fire period
## 94                inside - Complete data throughout fire period
## 95               outside - Complete data throughout fire period
## 97               outside - Complete data throughout fire period
## 99        outside - Sensor offline during fire, returned online
## 101              outside - Complete data throughout fire period
## 103              outside - Complete data throughout fire period
## 106        inside - Sensor offline during fire, returned online
## 115        inside - Sensor offline during fire, returned online
## 200  inside - Sensor offline during fire, did not return online
# map of damaged/destroyed buildings and AQ sensors (indoor/outdoor)
# chose not to do icons for destroyed/damaged buildings because the map got really busy
(fire_destruction_and_AQ_plot = leaflet(sensors) %>%
    addTiles() %>% 
    addAwesomeMarkers(icon = ~sensorIcon[Location],
                      lng = sensors$Lon, lat = sensors$Lat) %>% 
    addCircleMarkers(color = ~pal(destroyed_damaged$type),
                     radius = 3.5,
                     opacity = 1,
                     lng = destroyed_damaged$long, lat = destroyed_damaged$lat) %>%
    addLegend("topright", pal = pal, values = ~destroyed_damaged$type,
              title = "Fire damage"))
# create palette for location + time period
timeSensorIcon = awesomeIconList(
  "inside - Complete data throughout fire period" = makeAwesomeIcon(
    icon = "home",
    iconColor = "white",
    markerColor = "green",
    library = "fa"
  ),
  "inside - Sensor offline during fire, did not return online" = makeAwesomeIcon(
    icon = "tree",
    iconColor = "white",
    markerColor = "gray",
    library = "fa"
  ),
  "inside - Sensor offline during fire, returned online" = makeAwesomeIcon(
    icon = "home",
    iconColor = "white",
    markerColor = "blue",
    library = "fa"
  ),
   "outside - Complete data throughout fire period" = makeAwesomeIcon(
    icon = "home",
    iconColor = "white",
    markerColor = "green",
    library = "fa"
  ),
  "outside - Sensor offline during fire, did not return online" = makeAwesomeIcon(
    icon = "tree",
    iconColor = "white",
    markerColor = "gray",
    library = "fa"
  ),
  "outside - Sensor offline during fire, returned online" = makeAwesomeIcon(
    icon = "tree",
    iconColor = "white",
    markerColor = "blue",
    library = "fa"
  ))

# map based on when sensor was added
(fire_destruction_and_AQ_plot = leaflet(sensors) %>%
    addTiles() %>% 
    addAwesomeMarkers(icon = ~timeSensorIcon[loc_status],
                      lng = sensors$Lon, lat = sensors$Lat) %>% 
    addCircleMarkers(color = ~pal(destroyed_damaged$type),
                     radius = 3.5,
                     opacity = 1,
                     lng = destroyed_damaged$long, lat = destroyed_damaged$lat) %>%
    addLegend("topright", pal = pal, values = ~destroyed_damaged$type,
              title = "Fire damage") %>%
    addLegend("bottomright",
              pal = colorFactor(c("green","gray","blue"), domain = unique(sensors$Status)),
              values = ~sensors$Status, title = "Sensor Time Info"))

Create spatial objects

## Create spacetime objects -- 10 minute
# first, re-format datetime as a xts object
AQ_df$datetime = as.POSIXct(AQ_df$datetime)

# create spatial points object for each sensor
PA_sensors = SpatialPoints(AQ_df[!duplicated(AQ_df$ID), c("Lon", "Lat")],proj4string = CRS(prg))
summary(PA_sensors)
## Object of class SpatialPoints
## Coordinates:
##            min        max
## Lon -105.35808 -105.02234
## Lat   39.82646   40.13818
## Is projected: FALSE 
## proj4string : [+proj=longlat +datum=WGS84 +no_defs]
## Number of points: 144
# contruct spatiotemporal object
PA_STFDF = stConstruct(AQ_df, space = c("Lon","Lat"), time = "datetime", crs = CRS(prg), SpatialObj = sensors)

# turn ST object into a STFDF (stationary points)
PA_STFDF = as(PA_STFDF,"STFDF")

# see summary of data
summary(PA_STFDF)
## Object of class STFDF
##  with Dimensions (s, t, attr): (144, 17706, 12)
## [[Spatial:]]
## Object of class SpatialPoints
## Coordinates:
##            min        max
## Lon -105.35808 -105.02234
## Lat   39.82646   40.13818
## Is projected: FALSE 
## proj4string : [+proj=longlat +datum=WGS84 +no_defs]
## Number of points: 144
## [[Temporal:]]
##      Index                       timeIndex    
##  Min.   :2021-12-30 00:00:00   Min.   :    1  
##  1st Qu.:2022-01-29 17:42:30   1st Qu.: 4427  
##  Median :2022-03-01 11:25:00   Median : 8854  
##  Mean   :2022-03-01 11:25:00   Mean   : 8854  
##  3rd Qu.:2022-04-01 06:07:30   3rd Qu.:13280  
##  Max.   :2022-05-01 23:50:00   Max.   :17706  
## [[Data attributes:]]
##        X                 ID             Name             Location        
##  Min.   :      1   Min.   :  3062   Length:2549664     Length:2549664    
##  1st Qu.: 482877   1st Qu.: 82879   Class :character   Class :character  
##  Median : 965752   Median :110042   Mode  :character   Mode  :character  
##  Mean   : 965752   Mean   : 99993                                        
##  3rd Qu.:1448628   3rd Qu.:128801                                        
##  Max.   :1931503   Max.   :146952                                        
##  NA's   :618161    NA's   :618161                                        
##      pm25_a           pm25_b            temp              rh        
##  Min.   :   0.0   Min.   :   0.0   Min.   :-234.0   Min.   : 0.0    
##  1st Qu.:   0.0   1st Qu.:   0.0   1st Qu.:  38.2   1st Qu.:19.0    
##  Median :   0.5   Median :   0.4   Median :  50.0   Median :27.4    
##  Mean   :  34.6   Mean   :   6.6   Mean   :  50.2   Mean   :30.0    
##  3rd Qu.:   3.9   3rd Qu.:   3.7   3rd Qu.:  62.6   3rd Qu.:40.4    
##  Max.   :4999.0   Max.   :3542.5   Max.   : 127.4   Max.   :94.6    
##  NA's   :618161   NA's   :618161   NA's   :652435   NA's   :652435  
##   a_b_agree             pm_avg       ab_difference      corrected_pm   
##  Length:2549664     Min.   :   0.0   Min.   :-3537.2   Min.   :  -2.0  
##  Class :character   1st Qu.:   0.0   1st Qu.:   -0.1   1st Qu.:   3.3  
##  Mode  :character   Median :   0.6   Median :    0.0   Median :   4.2  
##                     Mean   :  20.6   Mean   :   28.0   Mean   :  18.6  
##                     3rd Qu.:   4.0   3rd Qu.:    0.2   3rd Qu.:   5.4  
##                     Max.   :2514.5   Max.   : 4994.6   Max.   :1978.1  
##                     NA's   :618161   NA's   :618161    NA's   :652435
# object class
class(PA_STFDF)
## [1] "STFDF"
## attr(,"package")
## [1] "spacetime"
# object dimensions
dim(PA_STFDF)
##     space      time variables 
##       144     17706        12

Map of all sensors we have data for (many not online during the)

ggplot(all_bounds, aes(geometry=geometry)) +
  geom_sf() +
  geom_sf(data=marshall_fire$geometry, fill="red", alpha=0.6) +
  geom_sf(data=st_as_sf(PA_sensors))

Time series plots

## Time series plots
stplot(PA_STFDF[,"2021-12-30::2021-12-31","temp"],mode="ts")

stplot(PA_STFDF[,"2021-12-30::2021-12-31","rh"],mode="ts")

stplot(PA_STFDF[,"2021-12-30::2021-12-31","corrected_pm"],mode="ts")

stplot(PA_STFDF[,"2021-12-30::2021-12-31","temp"],mode="xt")

stplot(PA_STFDF[,"2021-12-30::2021-12-31","rh"],mode="xt")

stplot(PA_STFDF[,"2021-12-30::2021-12-31","corrected_pm"],mode="xt")

Sensors With Complete Data

# make an STFDF with only the sensors that had complete data for the fire
complete_sensors = sensors %>% filter(Status == "Complete data throughout fire period")
complete_STFDF = PA_STFDF[PA_STFDF@data$ID %in% complete_sensors$ID]
stplot(complete_STFDF[,,"corrected_pm"], mode="ts")

stplot(complete_STFDF[,"2021-12-30::2021-12-31","corrected_pm"], mode="ts")

stplot(complete_STFDF[,"2021-12-30::2021-12-31","corrected_pm"], mode="tp")

Kriging

Funky Data from Variograms

# Index 15 -> ID 60517, STANDLEY LAKE
plot(fortify(PA_STFDF[PA_STFDF@data$ID == 60517]["2021-12-30::2022-01-02", "corrected_pm"]))

# looks normal

# Index 23 -> ID 91877, Eisenhower Dr
plot(fortify(PA_STFDF[PA_STFDF@data$ID == 91877]["2021-12-30::2022-01-02", "corrected_pm"]))

# has incredibly large values midday Friday - Sunday, could be snow related?

# Index 36 -> ID 112606, 1333 Wildwood Ct., Boulder, CO 80305 sensor
plot(fortify(PA_STFDF[PA_STFDF@data$ID == 112606]["2021-12-30::2022-01-02", "corrected_pm"]))

# looking at purpleair map, A channel is 'downgraded', 0% confidence

# Index 40 -> ID 119547, Broadlands sensor
plot(fortify(PA_STFDF[PA_STFDF@data$ID == 119547]["2021-12-30::2022-01-02", "corrected_pm"]))

# reasonable values

Starting with Dec 30 - Jan 2, only outside sensors.

# get all of the data for the 12/30 to 1/2 timeframe
krig_data <- complete_STFDF[,"2021-12-30::2022-01-02"]
# limit to only outdoor sensors
krig_data <- krig_data[krig_data@data$Location == "outside"]
# drop sensor 112606 b/c 0% confidence according to PurpleAir map
krig_data <- krig_data[krig_data@data$ID != 112606]
# drop values for "Eisenhower Dr" sensor that're over 1000
krig_data@data[(krig_data@data$ID == 91877) &
                 (as.numeric(krig_data@data$corrected_pm) > 1000) &
                 !is.na(krig_data@data$corrected_pm),]$corrected_pm <- NA
# get rid of the weird sensor that might be throwing residuals off
krig_data <- krig_data[krig_data@data$Name != "ntsky"]

# reproject the data to utm
reproj <- "+proj=utm +zone=13 +datum=WGS84 +units=km"
krig_data@sp <- spTransform(krig_data@sp, reproj)

# aggregate from 10 min intervals to hour intervals
krig_data <- aggregate(krig_data, by="hour", mean)

# split the data into the days we want
dec30 <- krig_data[, "2021-12-30"]
dec31 <- krig_data[, "2021-12-31"]
jan1 <- krig_data[, "2022-01-01"]
jan2 <- krig_data[, "2022-01-02"]

# split each day into the time periods of focus
# pre-fire period: betweeen 00:00 - 11:00 MST on 12/30/2021
pre_fire <- aggregate(dec30, by="11 hours", mean)[, "2021-12-30 00:00:00"]
# fire period (11am - 5pm)
during_fire <- aggregate(dec30[, "2021-12-30 11:00:00::2021-12-30 16:00:00"], by="6 hours", mean)
# post-fire on 12/30/2021
evening_of_fire <- aggregate(dec30[, "2021-12-30 16:00:00::2021-12-31 00:00:00"], by="8 hours", mean)
# dec 31st
december31 <- aggregate(dec31, by="day", mean)
# january 1st
january1 <- aggregate(jan1, by="day", mean)
# january 2nd
january2 <- aggregate(jan2, by="day", mean)

periods <- c(pre_fire, during_fire, evening_of_fire, december31, january1, january2)

# create variograms for these time periods
pre_fire.vgm <- variogram(corrected_pm~1, pre_fire[!is.na(pre_fire$corrected_pm),])
during_fire.vgm <- variogram(corrected_pm~1, during_fire[!is.na(during_fire$corrected_pm),])
evening_of_fire.vgm <- variogram(corrected_pm~1, evening_of_fire[!is.na(evening_of_fire$corrected_pm),])
december31.vgm <- variogram(corrected_pm~1, december31[!is.na(december31$corrected_pm),])
january1.vgm <- variogram(corrected_pm~1, january1[!is.na(january1$corrected_pm),])
january2.vgm <- variogram(corrected_pm~1, january2[!is.na(january2$corrected_pm),])

# plot the variograms
plot(pre_fire.vgm)

plot(during_fire.vgm)

plot(evening_of_fire.vgm)

plot(december31.vgm)

plot(january1.vgm)

plot(january2.vgm)

Variogram Clouds

for (period in periods) {
  cloud <- variogram(corrected_pm~1, period[!is.na(period$corrected_pm), ], cloud=TRUE)
  map <- variogram(corrected_pm~1, period[!is.na(period$corrected_pm), ], map=TRUE, cutoff=20, width=2)
  print(plot(cloud))
  print(plot(map))
}

No good models variogram models!!!!

pre_fire.fitted <- fit.variogram(pre_fire.vgm, vgm(c("Exp", "Sph", "Gau", "Mat")))
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit
## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit
## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit
pre_fire.fitted
##   model       psill    range
## 1   Nug    86.24564      0.0
## 2   Exp 62387.14334 -10613.7
plot(pre_fire.vgm, model=pre_fire.fitted)

##   model       psill    range
## 1   Nug    86.24564      0.0
## 2   Exp 62387.14334 -10613.7
during_fire.fitted <- fit.variogram(during_fire.vgm, vgm(c("Exp", "Sph", "Gau", "Mat")))
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?

## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?

## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
during_fire.fitted
##   model    psill    range
## 1   Nug 152.5208  0.00000
## 2   Exp 938.8172 68.95887
plot(during_fire.vgm, model=during_fire.fitted)

evening_of_fire.fitted <- fit.variogram(evening_of_fire.vgm, vgm(c("Exp", "Sph", "Gau", "Mat")))
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?

## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
evening_of_fire.fitted
##   model    psill    range
## 1   Nug   0.0000  0.00000
## 2   Gau 228.5243 23.99524
plot(evening_of_fire.vgm, model=evening_of_fire.fitted)

december31.fitted <- fit.variogram(december31.vgm, vgm(c("Exp", "Sph", "Gau", "Mat")))
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?

## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?

## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
december31.fitted
##   model     psill    range
## 1   Nug  1.141338  0.00000
## 2   Gau 38.915184 16.75872
plot(december31.vgm, model=december31.fitted)

january1.fitted <- fit.variogram(january1.vgm, vgm(c("Exp", "Sph", "Gau", "Mat")))
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
january1.fitted
##   model    psill    range
## 1   Nug 4.897684 0.000000
## 2   Sph 7.572974 2.524618
plot(january1.vgm, model=january1.fitted)

january2.fitted <- fit.variogram(january2.vgm, vgm(c("Exp", "Sph", "Gau", "Mat")))
january2.fitted
##   model     psill    range
## 1   Nug  1.572919  0.00000
## 2   Exp 33.290787 12.57439
plot(january1.vgm, model=january2.fitted)

Breaking down into smaller timeframes

for (i in 0:11) {
  if (i < 10) {
    i <- paste("0", i, sep="")
  }
  t <- paste("2021-12-30 ", i, ":00:00", sep="")
  i.data <- dec30[, t]
  i.vgm <- variogram(corrected_pm~1, i.data[!is.na(i.data$corrected_pm),])
  # print(plot(i.vgm, main=paste("Variogram at", t)))
  i.fitted <- fit.variogram(i.vgm, vgm(c("Exp", "Sph", "Gau", "Mat")))
  print(i.fitted)
  print(plot(i.vgm, model=i.fitted, main=paste("Variogram at", t)))
}
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
##   model    psill    range
## 1   Nug 0.000000 0.000000
## 2   Sph 3.904079 5.151113
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?

##   model    psill    range
## 1   Nug 0.000000 0.000000
## 2   Sph 3.095985 3.763504
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit
## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?

##   model    psill    range
## 1   Nug 0.000000 0.000000
## 2   Gau 3.260366 1.309588
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?

##   model     psill    range
## 1   Nug 0.4409823 0.000000
## 2   Sph 1.0171682 3.240835
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?

## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?

## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?

##   model     psill    range
## 1   Nug 0.1678131  0.00000
## 2   Gau 0.2917075 10.77585
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?

## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?

## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?

##   model     psill    range
## 1   Nug 0.0501557  0.00000
## 2   Gau 1.7184797 71.28375
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?

## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?

## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?

##   model    psill    range
## 1   Nug 0.028565  0.00000
## 2   Gau 0.201726 17.88329
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : linear model has singular covariance matrix

## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : linear model has singular covariance matrix
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?

## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?

##   model      psill   range
## 1   Nug 0.01918062 0.00000
## 2   Sph 0.01776751 1.34037
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?

## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?

## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?

## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?

##   model       psill    range
## 1   Nug 0.034425382  0.00000
## 2   Sph 0.005219032 29.42497
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit
## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit

##   model    psill     range
## 1   Nug  6181.04    0.0000
## 2   Exp 50155.36 -305.7408
##   model    psill     range
## 1   Nug  6181.04    0.0000
## 2   Exp 50155.36 -305.7408
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit
## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit

##   model    psill     range
## 1   Nug 391.9872    0.0000
## 2   Exp 375.1265 -189.0302
##   model    psill     range
## 1   Nug 391.9872    0.0000
## 2   Exp 375.1265 -189.0302
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit
## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?

##   model    psill    range
## 1   Nug 1.259568 0.000000
## 2   Gau 0.000000 4.613097

Still no good variograms!

Making grid for new data

grd <- SpatialPixels(SpatialPoints(as_Spatial(st_make_grid(all_bounds, n=50))), proj4string = prg)
plot(grd)

grd <- grd[as_Spatial(all_bounds),]
plot(grd)

grd <- spTransform(grd, CRS(reproj))

Plotting Corrected PM

This is just used to check the data from the sensors, and visualize if any of the sensors have odd values.

plot(pre_fire$corrected_pm)

plot(during_fire$corrected_pm)

plot(evening_of_fire$corrected_pm)

plot(december31$corrected_pm)

plot(january1$corrected_pm)

plot(january2$corrected_pm)

IDW

pre_fire.idw <- idw(corrected_pm~1, pre_fire[!is.na(pre_fire$corrected_pm),], grd)
## [inverse distance weighted interpolation]
spplot(pre_fire.idw[,1], main="Pre-Fire Period on 12/30")

during_fire.idw <- idw(corrected_pm~1, during_fire[!is.na(during_fire$corrected_pm),], grd)
## [inverse distance weighted interpolation]
spplot(during_fire.idw[,1], main="During Fire on 12/30")

evening_of_fire.idw <- idw(corrected_pm~1, evening_of_fire[!is.na(evening_of_fire$corrected_pm),], grd)
## [inverse distance weighted interpolation]
spplot(evening_of_fire.idw[,1], main="After Fire on 12/30")

december31.idw <- idw(corrected_pm~1, december31[!is.na(december31$corrected_pm),], grd)
## [inverse distance weighted interpolation]
spplot(december31.idw[,1], main="December 31st")

january1.idw <- idw(corrected_pm~1, january1[!is.na(january1$corrected_pm),], grd)
## [inverse distance weighted interpolation]
spplot(january1.idw[,1], main="January 1st")

january2.idw <- idw(corrected_pm~1, january2[!is.na(january2$corrected_pm),], grd)
## [inverse distance weighted interpolation]
spplot(january2.idw[,1], main="January 2nd")

aqi_colors <- c("#00e400",
                "#ffff00",
                "#ff7e00",
                "#ff0000",
                "#8f3f97",
                "#7e0023")

aqi_labels <- c("Good (0-12.0)",
                "Moderate (12.1-35.4)",
                "Unhealthy for Sensitive Groups (35.5-55.4)",
                "Unhealthy (55.5-150.4)",
                "Very Unhealthy (150.5-250.4)",
                "Hazardous (250.5-500.4)")

factorize_aqi <- function (idw) {
  idw$PM2.5 <- cut(idw$var1.pred, breaks=c(0, 12.0, 35.4, 55.4, 150.4, 250.4), na.omit=T)
  return(idw)
}

plot_idw <- function(idw, title) {
  return(ggplot(all_bounds) +
  geom_sf() +
  geom_sf(data=st_as_sf(factorize_aqi(idw)), aes(col=PM2.5), size=2) +
  scale_color_manual(values=aqi_colors, labels=aqi_labels, name="PM2.5 (µg/m^3)") +
  geom_sf(data=marshall_fire$geometry, col="gray20", alpha=0.1) +
  ggtitle(title))
}

plot_idw(pre_fire.idw, "Pre Fire IDW")

plot_idw(during_fire.idw, "During Fire IDW")

plot_idw(evening_of_fire.idw, "Evening of fire IDW")

plot_idw(december31.idw, "December 31st IDW")

plot_idw(january1.idw, "January 1st IDW")

plot_idw(january2.idw, "January 2nd IDW")

r_grd <- raster(grd, nrows=50, ncols=49)

# factorize_aqi <- function (idw) {
#   idw$PM2.5 <- cut(idw$var1.pred, breaks=c(0, 12.0, 35.4, 55.4, 150.4, 250.4), na.omit=T)
#   return(idw)
# }

# scale_color_manual(values=aqi_colors, labels=aqi_labels, name="PM2.5 (µg/m^3)") +

plot_idw_raster <- function(idw, title) {
  r <- raster::rasterize(idw, r_grd, field="var1.pred")
  extent(r) <- st_bbox(all_bounds)
  r_df <- r %>%
    rasterToPoints() %>%
    data.frame() %>%
    mutate(cuts = cut(layer, breaks=c(0, 12.0, 35.4, 55.4, 150.4, 250.4), na.omit=T))
  
  return(ggplot() +
          geom_tile(data=r_df, aes(x=x, y=y, fill=cuts))+
          scale_fill_manual(values=aqi_colors, labels=aqi_labels, name="PM2.5 (µg/m^3)") +
          geom_sf(data=marshall_fire$geometry, col="gray20", alpha=0.1) +
          ggtitle(title))
}

plot_idw_raster(pre_fire.idw, "Pre Fire IDW")

plot_idw_raster(during_fire.idw, "During Fire IDW")

plot_idw_raster(evening_of_fire.idw, "Evening of fire IDW")

plot_idw_raster(december31.idw, "December 31st IDW")

plot_idw_raster(january1.idw, "January 1st IDW")

plot_idw_raster(january2.idw, "January 2nd IDW")

Cross Validation

pre_fire.idw_cv <- krige.cv(corrected_pm~1, pre_fire[!is.na(pre_fire$corrected_pm),], nfold=38)
pre_fire.idw_sd <- sd(pre_fire.idw_cv$residual)
paste("Pre fire standard deviation:", pre_fire.idw_sd)
## [1] "Pre fire standard deviation: 7.58649993004707"
spplot(pre_fire.idw_cv, main="Pre-Fire Period on 12/30")

pre_fire.idw_cv$zscore <- pre_fire.idw_cv$residual / pre_fire.idw_sd

during_fire.idw_cv <- krige.cv(corrected_pm~1, during_fire[!is.na(during_fire$corrected_pm),], nfold=38)
during_fire.idw_sd <- sd(during_fire.idw_cv$residual)
paste("During fire standard deviation:", during_fire.idw_sd)
## [1] "During fire standard deviation: 15.376531044756"
spplot(during_fire.idw_cv, main="During Fire on 12/30")

during_fire.idw_cv$zscore <- during_fire.idw_cv$residual / during_fire.idw_sd

evening_of_fire.idw_cv <- krige.cv(corrected_pm~1, evening_of_fire[!is.na(evening_of_fire$corrected_pm),], nfold=38)
evening_of_fire.idw_sd <- sd(evening_of_fire.idw_cv$residual)
paste("Evening of fire standard deviation:", evening_of_fire.idw_sd)
## [1] "Evening of fire standard deviation: 4.85088863186492"
spplot(evening_of_fire.idw_cv, main="After Fire on 12/30")

evening_of_fire.idw_cv$zscore <- evening_of_fire.idw_cv$residual / evening_of_fire.idw_sd

december31.idw_cv <- krige.cv(corrected_pm~1, december31[!is.na(december31$corrected_pm),], nfold=38)
december31.idw_sd <- sd(december31.idw_cv$residual)
paste("December 31st standard deviation:", december31.idw_sd)
## [1] "December 31st standard deviation: 2.36672112753216"
spplot(december31.idw_cv, main="December 31st")

december31.idw_cv$zscore <- december31.idw_cv$residual / december31.idw_sd

january1.idw_cv <- krige.cv(corrected_pm~1, january1[!is.na(january1$corrected_pm),], nfold=38)
january1.idw_sd <- sd(january1.idw_cv$residual)
paste("January 1 standard deviation:", january1.idw_sd)
## [1] "January 1 standard deviation: 3.68472999489718"
spplot(january1.idw_cv, main="January 1st")

january1.idw_cv$zscore <- january1.idw_cv$residual / january1.idw_sd

january2.idw_cv <- krige.cv(corrected_pm~1, january2[!is.na(january2$corrected_pm),], nfold=38)
january2.idw_sd <- sd(january2.idw_cv$residual)
paste("January 2 standard deviation:", january2.idw_sd)
## [1] "January 2 standard deviation: 3.6116917414784"
spplot(january1.idw_cv, main="January 2nd")

january2.idw_cv$zscore <- january2.idw_cv$residual / january2.idw_sd
plot_idw_resid <- function(idw.cv, title) {
  return(ggplot(all_bounds) +
  geom_sf(fill="#FbFbFb", color="gray") + 
  geom_sf(data=marshall_fire$geometry, col="red", alpha=0.1) +
  geom_sf(data=st_as_sf(idw.cv), aes(fill=residual), color="black", shape=21, size=2) +
  scale_color_gradient2(aesthetics="fill") +
  ggtitle(title))
}

plot_idw_resid_zscore <- function(idw.cv, title) {
  return(ggplot(all_bounds) +
  geom_sf(fill="#FbFbFb", color="gray") + 
  geom_sf(data=marshall_fire$geometry, col="red", alpha=0.1) +
  geom_sf(data=st_as_sf(idw.cv), aes(fill=zscore), color="black", shape=21, size=2) +
  scale_color_gradient2(aesthetics="fill") +
  ggtitle(title))
}

plot_idw_cv <- function(idw, idw_title, idw.cv) {
  p1 <- plot_idw(idw, idw_title)
  p2 <- plot_idw_resid(idw.cv, "Residuals")
  return(grid.arrange(p1, p2, ncol=2))
}

plot_idw_resid(pre_fire.idw_cv, "Pre-Fire Residuals")

plot_idw_resid(during_fire.idw_cv, "During Fire Residuals")

plot_idw_resid(evening_of_fire.idw_cv, "After Fire Residuals")

plot_idw_resid(december31.idw_cv, "December 31st Residuals")

plot_idw_resid(january1.idw_cv, "January 1st Residuals")

plot_idw_resid(january2.idw_cv, "January 2nd Residuals")

plot_idw_resid_zscore(pre_fire.idw_cv, "Pre-Fire Residual Z-Scores")

plot_idw_resid_zscore(during_fire.idw_cv, "During Fire Residual Z-Scores")

plot_idw_resid_zscore(evening_of_fire.idw_cv, "After Fire Residual Z-Scores")

plot_idw_resid_zscore(december31.idw_cv, "December 31st Residual Z-Scores")

plot_idw_resid_zscore(january1.idw_cv, "January 1st Residual Z-Scores")

plot_idw_resid_zscore(january2.idw_cv, "January 2nd Residual Z-Scores")

plot_idw_cv(pre_fire.idw, "Pre-Fire IDW", pre_fire.idw_cv)

plot_idw_cv(during_fire.idw, "During Fire IDW", during_fire.idw_cv)

plot_idw_cv(evening_of_fire.idw, "Evening of Fire IDW", evening_of_fire.idw_cv)

plot_idw_cv(december31.idw, "December 31st IDW", december31.idw_cv)

plot_idw_cv(january1.idw, "January 1st IDW", january1.idw_cv)

plot_idw_cv(january2.idw, "January 2nd IDW", january2.idw_cv)

Notes:

  • Pre-fire period the residuals look fairly good, some over & underprediction in downtown Boulder
  • During fire period has 1 large underprediction value, few other overpredicted near fire area
  • Evening after the fire has underprediction, but not much
  • December 31st has high under & over predicted values right next to each other near fire
  • Same as December 31st on Jan 1st
  • Same as prior 2 days
    • Something happened starting on December 31st with these two sensors? Will look into both and see how the values differ over time
# One sensor is "Eisenhower Dr"
plot(fortify(PA_STFDF[PA_STFDF@data$Name == "Eisenhower Dr"]["2021-12-30::2022-01-02", "corrected_pm"]))

# The other is "Purple House"
plot(fortify(PA_STFDF[PA_STFDF@data$Name == "Purple House"]["2021-12-30::2022-01-02", "corrected_pm"]), col="red")

The issue is coming from the Eisenhower Dr sensor most likely, as this is the one that seems messed up by the snow. Going to see how everything looks excluding the high values for this sensor.

plot(fortify(PA_STFDF[PA_STFDF@data$Name == "Eisenhower Dr"]["2021-12-30::2022-01-02", "corrected_pm"]) %>%
               filter(corrected_pm <= 1000))

Checking the sensor on the edge of the fire boundary with high residuals/z scores.

plot(fortify(PA_STFDF[PA_STFDF@data$Name == "ntsky"]["2021-12-30::2022-01-02", "corrected_pm"]))

Looping at Hourly Timescale

for (j in c(30, 31, 01, 02)) {
  if (j > 3) {
    date <- paste("2021-12-", j, sep="")
  } else {
    date <- paste("2022-01-", j, sep="")
  }
  for (i in 0:23) {
    if (i < 10) {
      i <- paste("0", i, sep="")
    }
    t <- paste(date, " ", i, ":00:00", sep="")
    i.data <- krig_data[, t]
    i.idw <- idw(corrected_pm~1, i.data[!is.na(i.data$corrected_pm),], grd)
    print(plot_idw_raster(i.idw, paste(t, "IDW")))
  }
}

for (j in c(30, 31, 01, 02)) {
  if (j > 3) {
    date <- paste("2021-12-", j, sep="")
  } else {
    date <- paste("2022-01-", j, sep="")
  }
  for (i in 0:23) {
    if (i < 10) {
      i <- paste("0", i, sep="")
    }
    t <- paste(date, " ", i, ":00:00", sep="")
    i.data <- krig_data[, t]
    i.idw <- idw(corrected_pm~1, i.data[!is.na(i.data$corrected_pm),], grd)
    i.idw_cv <- krige.cv(corrected_pm~1, i.data[!is.na(i.data$corrected_pm),], nfold=38)
    i.idw_sd <- sd(i.idw_cv$residual)
    i.idw_cv$zscore <- i.idw_cv$residual / i.idw_sd
    print(plot_idw_resid_zscore(i.idw_cv, paste(t, "IDW Residual Z-Scores")))
  }
}

Daily Gifs Starting 1/3/2022

# get the full timerange of data & clean it up
daily <- complete_STFDF[,"2022-01-01::"]
daily <- daily[daily@data$Location == "outside"]
# drop sensor 112606 b/c 0% confidence according to PurpleAir
daily <- daily[daily@data$ID != 112606]
# drop values for "Eisenhower Dr" sensor that're over 1000
daily@data[(daily@data$ID == 91877) &
                 (as.numeric(daily@data$corrected_pm) > 1000) &
                 !is.na(daily@data$corrected_pm),]$corrected_pm <- NA
# get rid of the weird sensor that might be throwing residuals off
daily <- daily[daily@data$Name != "ntsky"]
# get rid of "ponderosa" sensor values after april 11 (when sensor looks to get wonky)
daily@data[(daily@data$ID == 102614) & (daily@endTime < "2022-04-11") & !is.na(daily@data)] <- NA

# reproject & aggregate
daily@sp <- spTransform(daily@sp, reproj)
daily <- aggregate(daily, by="day", mean)
for (j in 1:4) {
  date <- paste("2022-0", j, sep="")
  if (j %in% c(1, 3, 5)) {
    range = 1:31
  } else if (j == 2) {
    range = 1:28
  } else {
    range = 1:30
  }
  
  for (i in range) {
    if (i < 10) {
      i <- paste("0", i, sep="")
    }
    day <- paste(date, "-", i, sep="")
    
    data <- daily[, day]
    idw <- idw(corrected_pm~1, data[!is.na(data$corrected_pm),], grd)
    print(plot_idw_raster(idw, paste(day, "IDW")))
    }
}

plot(daily[, "2022-04-25"]$corrected_pm) # sensor at index 29 is weird = sensor ID 102614

plot(fortify(PA_STFDF[PA_STFDF@data$ID == 102614]['2022-04', "corrected_pm"]))

# ponderosa sensor has 0% confidence on purpleair map at present date
# png(file = "images/month_map.png", width = 6000, height = 3800, res = 300)
# (month_plot = basemap +
#   geom_point(data = month_added_plot, aes(x=Lon,y=Lat, label = ID, shape = Location, color = Month)) +
#    facet_grid(cols = vars(Month)))
# dev.off()